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■ We present a new algorithm to sample the constrained eigenvalues of the ini- 

tial shear field associated with Gaussian statistics, called the 'peak/dip excursion-set- 
based' algorithm, at positions which correspond to peaks or dips of the correlated 
density field. The computational procedure is based on a new formula which extends 
Doroshkevich's unconditional distribution for the eigenvalues of the linear tidal field, 
to account for the fact that halos and voids may correspond to maxima or minima 
of the density field. The ability to differentiate between random positions and special 
points in space around which halos or voids may form (i.e. peaks/dips), encoded in 
the new formula and reflected in the algorithm, naturally leads to a straightforward 
implementation of an excursion set model for peaks and dips in Gaussian random 

■ fields - one of the key advantages of this sampling procedure. In addition, it offers 
' novel insights into the statistical description of the cosmic web. As a first physical ap- 
, plication, we show how the standard distributions of shear ellipticity and prolatcncss 

in triaxial models of structure formation arc modified by the constraint. In particular, 
we provide a new expression for the conditional distribution of shape parameters given 
, the density peak constraint, which generalizes some previous literature work. The for- 

mula has important implications for the modeling of non-spherical dark matter halo 
shapes, in relation to their initial shape distribution. We also test and confirm our 
theoretical predictions for the individual distributions of eigenvalues subjected to the 
k> ' extremum constraint, along with other directly related conditional probabilities. Fi- 

j_j \ nally, we indicate how the proposed sampling procedure naturally integrates into the 

C$ ■ standard excursion set model, potentially solving some of its well-known problems, 

and into the ellipsoidal collapse framework. Several other ongoing applications and 
extensions, towards the development of algorithms for the morphology and topology 
of the cosmic web, are discussed at the end. 

Key words: methods: analytical - methods: statistical - galaxies: formation - galax- 
ies: halos - cosmology: theory - large-scale structure of Universe. 



1 INTRODUCTION 

The large-scale spatial organization of matter in clusters, filaments, sheets and voids, known as the cosmic web (Peebles 1980; 
Bardeen et al. 1986; Bond et al. 1991), is the manifestation of the anisotropic nature of gravitational collapse. This typical 
filamentary pattern has been confirmed by observations, for example with the 2dFGRS, SDSS and 2MASS redshift surveys of 
the nearby Universe (Colless et al. 2003; Tegmark et al. 2004; Huchra et al. 2005), and is routinely seen in large scale A-body 
numerical simulations - see for example the New Horizon Runs (Kim et al. 2011) or the recent Millennium-XXL (Angulo et 
al. 2012). Basic characteristics of the cosmic web are the spatial arrangement of galaxies and mass into elongated filaments, 
sheet-like walls and dense compact clusters, the existence of large near-empty void regions, and the hierarchical nature of this 
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mass distribution (Aragon-Calvo et al. 2007, 2010a,b; Zhang et al. 2009). In particular, as pointed out by Bond, Kofman & 
Pogosyan (1996), 'embryonic' cosmic web is already present in the primordial density field. These key properties, along with 
the alignment of shape and angular momentum of objects (i.e. Argyres et al. 1986; Catelan et al. 2001; Lee & Springel 2010), 
are mainly due to the effects of the tidal field, associated with the gravitational potential - while the Hessian of the density 
field (i.e. the inertia tensor) plays a secondary role in determining the characteristic pattern of the cosmic web. 

Pioneering works devoted to the key role of the initial tidal field in shaping large-scale structures trace back to Doroshke- 
vich & Zeldovich (1964), Doroshkevich (1970), Zeldovich (1970), Sunyaev & Zeldovich (1972), Icke (1973), and Doroshkevich 
& Shandarin (1978); their studies have contributed to reach a solid understanding of the formation and evolution of struc- 
tures. Other classical works in this direction (i.e. Peebles 1980; White 1984; Bardeen et al. 1986; Kaiser 1986; Bertschinger 
1987; Dubinski 1992; Bond & Myers 1996; Bond, Kofman & Pogosyan 1996; van de Weygaert & Bertschinger 1996) have 
considerably improved the statistical description of the cosmic web from first principles. Their impact is broader, as there 
is a correspondence between structures in the evolved density field and local properties of the linear tidal shear; this allows 
one to estimate the morphology of the cosmic web (Bond & Myers 1996; Rossi, Sheth & Tormen 2011), and is crucial in 
understanding the nonlinear evolution of cosmic structures (Springel et al. 2005; Shandarin et al. 2006; Pogosyan et al. 2009), 
the statistical properties of voids (Lee & Park 2006; Platen, van de Weygaert & Jones 2008; Lavaux & Wandelt 2010), the 
alignment of shape and angular momentum of halos (West 1989; Catelan et al. 2001; Faltenbacher et al. 2009; Lee & Springel 
2010), and for characterizing the geometry and topology of the cosmic web (Gott et al. 1986, 1989; Park & Gott 1991; Park 
et al. 1992, 2005; van de Weygaert & Bond 2008; Forero-Romero et al. 2009; Aragon-Calvo et al. 2010a,b; Choi et al. 2010; 
Shandarin et al. 2010; van de Weygaert et al. 2011; Cautun et al. 2012; Hidding et al. 2012). In addition, the eigenvalues of 
the mass tidal tensor can be used to classify the large-scale environment (Shen et al. 2006; Hahn et al. 2007; Zhang et al. 
2009) - hence their fundamental importance. 

In the context of the initial shear field, Doroshkevich (1970) provided a key contribution by deriving the joint probability 
distribution of an ordered set of eigenvalues in the tidal field matrix corresponding to a Gaussian potential, given the variance 
of the density field. Throughout the paper, we will refer to it as the unconditional probability distribution of shear eigenvalues. 
Because the initial shear field associated with Gaussian statistics plays a major role in the formation of large scale structures, 
considerable analytic work has been based on Doroshkevich's unconditional formulas since their appearance, but those relations 
neglect the fact that halos (voids) may correspond to maxima (minima) of the underlying density field. The local extrema 
of such field (i.e. peaks/dips) are plausible sites for the formation of nonlinear structures (Bardeen et al. 1986), and some 
numerical studies have indeed reported a good correspondence between peaks in the initial conditions and halos at late times 
(see in particular Ludlow & Porciani 2011). Hence, their statistical properties can be used to predict the abundances and 
clustering of objects of various types, and in studies of triaxial formation of large-scale structures (Bardeen et al. 1986; Bond 
& Myers 1996). 

Motivated by these reasons, recently Rossi (2012) has provided a set of analytic expressions which extend the work of 
Doroshkevich (1970) and Bardeen et al. (1986) - and are akin in philosophy to that of van de Weygaert & Bertschinger (1996). 
These new relations incorporate the peak/dip constraint into the statistical description of the initial shear field, and are able 
to differentiate between random positions and peak/dips in the correlated density field. In essence, they allow one to express 
the joint probability distribution of an ordered set of eigenvalues in the initial shear field, given the fact that positions are 
peaks or dips of the density - and not just random spatial locations. These relations are obtained by requiring the density 
Hessian (i.e. the matrix of the second derivatives of the density field, associated with the curvature) to be positive/negative 
definite, which is the case in the vicinity of minima/maxima of the density. The correlation strength between the gravitational 
and density fields is quantified via a reduced parameter 7, which plays a major role in peaks theory (i.e. 7 is the same as in 
Bardeen et al. 1986 for Gaussian smoothing filters - see Appendix [A). Doroshkevich's (1970) unconditional formulas are then 
naturally recovered in the absence of correlation, when 7 = 0. From these new conditional joint probabilities, it is possible 
to derive the individual distributions of eigenvalues subjected to the extremum constraint, along with some other related 
conditional probabilities; their expressions were also provided in Rossi (2012), extending the work of Lee & Shandarin (1998). 

The primary goal of this paper is show how the main conditional formulas presented in Rossi (2012) and reported here 
(Equations HI andl lip lead to a new algorithm - called the 'peak/dip excursion-set-based' algorithm - to sample the constrained 
eigenvalues of the initial shear field associated with Gaussian statistics, at positions which correspond to peaks or dips of the 
correlated density field. While it is clearly possible to sample the constrained eigenvalues of the tidal field directly from the 
conditional probability distribution function, as done for example by Lavaux & Wandelt (2010) in the context of cosmic voids, 
the theoretical work carried out by Rossi (2012) allows for a much simpler procedure, part of which was previously thought 
not to be achievable analytically (see again the Appendix B in Lavaux & Wandelt 2010). Besides providing novel theoretical 
insights, the main strength of the algorithm resides in its natural inclusion into the standard excursion set framework, allowing 
for a generalization. As we will discuss in Section [6] this technique potentially solves a long-standing problem of the standard 
excursion set theory; moreover, it is well-suited for the ellipsoidal collapse framework. 

The other goal of the paper is to present a first physical application of the new algorithm, related to the morphology of 
the cosmic web. In particular, we provide a new expression for the conditional distribution of shape parameters (i.e. ellipticity 
and prolateness) in the presence of the density peak constraint, which generalizes some previous literature work and combines 
the formalism of Bardeen et al. (1986) - based on the density field - with that of Bond & Myers (1996) - based on the shear 
field. The formula has important implications for the modeling of non-spherical dark matter halos and their evolved halo 
shapes, in relation to the initial shape distribution. In addition, we also test and confirm our theoretical predictions for the 
individual distributions of eigenvalues subjected to the extremum constraint, along with other directly related conditional 
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probabilities. Finally, we illustrate how this algorithm can be readily merged into the excursion set framework (Peacock & 
Heavens 1990; Bond et al. 1991), and in particular to obtain an excursion set model for peaks and dips in Gaussian random 
fields. The key point is the ability to differentiate between random positions and peaks/dips, which is contained in Equations 
|T} and (111[) and encapsulated in the algorithm. While we leave this latter part to a dedicated forthcoming publication, we 
anticipate the main ideas here. We also discuss several other ongoing applications and extensions, towards the development 
of algorithms for classifying cosmic web structures. 

The layout of the paper is organized as follows. Section [2] provides a short review of the key equations derived in 
Rossi (2012), which constitute the theoretical framework for the new algorithm described here; the main notation adopted 
is summarized in Appendix [X] for convenience. Moving from this mathematical background, Section [3] presents the new 
'peak/dip excursion-set-based' algorithm (some insights derived from this part are left in Appendix IB) , while Section [4] tests 
its performance against various analytic distributions of eigenvalues and related conditional probabilities, subjected to the 
extremum constraint. Section [S] shows a physical application of the computational procedure, towards the morphology of the 
cosmic web. A new expression for the conditional distribution of ellipticity and prolateness in the presence of the density 
peak constraint is also given; in particular, the description of Bardeen et al. (1986) is combined with that of Bond & Myers 
(1996). Section HJ] illustrates how this new algorithm can be readily used to implement an excursion set model for peaks and 
dips in Gaussian random fields, and makes the connection with some previous literature. Finally, Section [7] discusses several 
ongoing and future promising applications, which will be presented in forthcoming publications, and in particular the use of 
this algorithm for triaxial models of collapse and in relation to the morphology and topology of the cosmic web. 



2 JOINT CONDITIONAL DISTRIBUTION OF EIGENVALUES IN THE PEAK/DIP PICTURE 

We begin by reexamining two main results derived in Rossi (2012), which constitute the key for developing a new algorithm 
to sample the constrained eigenvalues of the initial shear field. This part may be also regarded as a compact review of the 
main formulas for the constrained shear eigenvalues, which can be used directly for several applications related to the cosmic 
web. The notation adopted here is the same as the one introduced by Rossi (2012), with a few minor changes to make the 
connection with previous literature more explicit. It is summarized in Appendix [X] a reader not familiar with the notation 
may want to start from the appendix first. In particular, in what follows we do not adopt 'reduced' variables, so that the 
various dependencies on a values (i.e. ctt = (To for the shear, and cth = ct 2 for the density Hessian) are now shown explicitly. 
However, for the sake of clarity, we omit to indicate the understood dependence of these two global parameters in the left-hand 
side of all the formulas. Note that we prefer to write ctt and cth, rather than cto and 02, for their more intuitive meaning 
(i.e. the label T indicates that a quantity is connected to the shear field, while the label H points to the Hessian of the 
density field). As explained in Appendix lAl we also introduce the 6-dimensional vectors T = (Tn, T22, T33, Ti 2 , T13, T23) and 
H = (Hn, H22, H33, H\2, H13, H23), derived from the components of their corresponding shear and density Hessian tensors. 

Rossi (2012) obtained the following expression for the probability of observing a tidal field T for the gravitational potential, 
given a curvature H for the density field: 



P(TiH <^ = 16^3CT|(l- 72 )3 ^[-^ T (W)(^ 2 " H « 

where 

Ki = (Tn - r)Hxi) + (T 22 - VH22) + (T 33 - VH33) = kx - r)hx (2) 

K 2 = (Tu -T?ffii)(T 22 - 77.ff 22 ) + (Tu -^11)^33-77^33) 

+ (T 22 - r)H 22 ){T 3 3 - rjHss) - (T12 - vHi 2 f - (T13 - vHisf - (T23 - VH23) 2 = k 2 + v 2 h 2 - rjhxkx + j?t (3) 

t = TiiHn + T22H22 + T33H33 + 2T12H12 + 2T13H13 + 2T23H23 (4) 

fci = Tn + T22 + T33 (5) 

k2 = T1T22 + T11T33 + T22T33 — Ty 2 — T13 — T 2 2 3 (6) 

hx = Hxx+H 22 + H 3 3 (7) 

h.2 = H11H22 + H11H33 + H22H33 — H12 — H13 — H23 (8) 

r\ = 7ctt/cth. (9) 



The corresponding unconditional marginal distributions p(T) and p(H) are multidimensional Gaussians, expressed using 
Doroshkevich's formulas as 

1 c3 -1 r o -, 1 c l 3 1r3 1 

P(T) = kVI^ 4 exp I^ (2k? - 5k2) l ' P(H) = I^^ exp [^ (2h? - 5h2) l • (10) 

Equation {TJ generalizes Doroshkevich's formulas (|10|) to include the fact that halos/voids may correspond to maxima/minima 
of the density field. Note in particular that p(T|H, 7) is a multivariate Gaussian distribution with mean b — rjil and covariance 
matrix (1 — 7 2 )ct t A/15, with A given in Appendix(A] Clearly, one can also consider the reverse distribution p(H|T,7), the 
expression of which is given in Rossi (2012). This is useful in order to make the connection and extend some results derived 
in Bardeen et al. (1986), the subject of a forthcoming publication. 
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It is also possible to express ([T]) in terms of the constrained eigenvalues of T|H, indicated with (j (i = 1, 2, 3) and ordered 
as Ci ^ C2 ^ C3- The result is: 



1^1 r 3 1 

p(Cl ' C2 ' C3|7) = iV5* - ^ ° XP [ " 2^(1 - 7 2 ) (2 ^ - 5 * a) l (Cl - Ca)(Cl - C3)(Ca " C3) 

where in terms of constrained eigenvalues we now have: 

#i = Ci +C2 + C3 =ki -rjhi (12) 

#2 = C1C2 + C1C3 + C2C3 = k 2 +r) 2 h 2 -rjhiki +7r (13) 

r = Ai£i + A2C2 + A 3 £ 3 (14) 

fci = Ai + A 2 + A 3 (15) 

k 2 = A1A2 + A1A3 + A2A3 (16) 

hi = £1+6 + 6 (17) 

fca = fik+fifa + kk (18) 

Ci = Ai-^i. (19) 

The partial distributions p(Ai, A2, A3) and p(£i, £2, £3) are expressed by Doroshkevich's unconditional formulas as 

p(Ai,A 2 ,A 3 ) = -ij-^exp[--^(2k? - 5k 2 )l(Ai - A 2 )(Ai - A 3 )(A 2 - A3) (20) 

8V 07T & T 1 J 



15 1 r 3 1 

P(£i , 6, 60 = -7=- -r exp - — 3-(2h? - 5h 2 ) (Ci - 6) (6 - 6) (6 - 6), (21) 

where Ai,A2,A3 are the eigenvalues of the shear tensor, while £1,^2, £3 are those of the density Hessian. Similarly, one can 
easily obtain formulas for the reverse probability functions. Note also that ki — Ai + A2 + A3 is simply the overdensity 5t 
associated to the shear field, while £1 + £2 + £3 = hi = <5h- Later on, we will make use of the peak height v = <5t/ot and of 
the peak curvature x = Sn/cra- 

Equations |l} and (|lip allow one to develop a new algorithm to sample the constrained eigenvalues of the initial shear 
field, presented in the next section. It will be then straightforward to use this algorithm in order to implement an excursion 
set model for peaks and dips in Gaussian random fields. 



3 THE PEAK/DIP EXCURSION-SET-BASED ALGORITHM 

While the constrained eigenvalues of the initial shear field can be sampled directly from their probability distribution function 
(i.e. Equations HI or [TT|l . the new theoretical formalism developed by Rossi (2012) leads to a much simpler algorithm, which 
is particularly interesting because well-suited for the ellipsoidal collapse model (Section [5]), and naturally integrable in the 
excursion set framework (Section|6]). In addition, the algebra leading to the new computational procedure allows one to explain 
analytically how the halo (void) shape distributions are altered by the inclusion of the peak (dip) constraint, and offers a 
variety of applications and insights on the statistical description of the cosmic web (see Section [7] and Appendix [B]). In what 
follows, we present the mathematical aspects of the 'peak/dip excursion-set-based' algorithm, and describe in detail the novel 
computational procedure. 

In particular, we are mainly interested in the distribution p(T|H > 0,7), the joint probability of observing a tidal field 
T for the gravitational potential with a positive density curvature H > (i.e. at density peak locations). Clearly, 

„ (T |H > 7) - PC^XW _ Jh>qP(H)p(T|H, 7 ) dH 

p(l|H>0, 7 )- p(H>Q) - /H>oP(H )dH ' [22) 

where p(T|H, 7) and p(H) are given by Equations {TJ and (|10[) . and the integrals are 6-dimensional. A similar expression can 
be written in terms of the constrained distributions of eigenvalues, using Equation (|lll) instead. In principle, one should then 
compute the previous integral to obtain p(T|H > 0,7). However, it is easier to sample p(T|H, 7 ) and impose the condition 
H > directly, so that we are effectively computing p(T|H > 0, 7). Given the nature of p(T|H,7) and p(H) expressed by {1} 
and pop , this can be readily achieved - as the following algebra will show. 

In fact, because the elements of the density Hessian are drawn from a multivariate Gaussian distribution with zero mean 
and covariance matrix <thA/15, where A is given by Equation (|A8[) . one can simply obtain them by generating six independent 
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zero- mean unit-variance Gaussian random variates yi (i = 1,6), and then determine the various components as: 





0"H ( 

in 






H22 


0"H / 
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H33 


0"H / 
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[- yi -±y2 + 
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H\2 


= H21 


OH 

= 7ff y4 




H13 


= H31 


= 7ff y5 




H23 


= H32 







(23) 

Similarly, the elements of the shear tensor are also drawn from a multivariate Gaussian distribution with mean zero and 
covariance matrix <ttA/15. Hence, if z\ (i — 1, 6) are other six independent zero-mean unit-variance Gaussian random variates, 
the components of the shear tensor are given by: 

T„ = + 

n. - =(-.- 

T12 = T21 = g4 







TT 2 ) 




1 


CO 


-^22 - 


"7r 3 


1 


CO 


75 


7= z 3 



15 

T13 = r 3 i = -^Lz 5 



15 

t 23 = r 32 = 4Lz 6 . (24) 

Next, consider the 6-dimensional vector T|H, made from the components of the conditional shear field. Its elements are 
drawn from a multivariate Gaussian distribution p(T|H,7) with mean b = r/H and covariance matrix (1 — j 2 )a^A/15, where 
the elements of the density Hessian are expressed by (|23p . Therefore, this implies that the components of T|H are distributed 
according to: 



T22\H22 = VH22 + (-i 1 -_ fa -_ fc j =T (- mi -_ ma -_ m8 



T 33 \H 33 = ^33 + ^ 3 72 (-/ 1 -^ 2 + ^3)=f (-rm-±m 2 + ^=m 3 



— 7 



2 



v/15 715 



15 V15 



— 7 



2 



T23IH23 = 7321^32 = ^23+ V ^ l 6 =^=m 6 . (25) 

V15 Vl5 



In the previous expression, the various Zj (i = 1, 6) are other six independent zero-mean unit-variance Gaussian variates, while 
the mi (i — 1,6) are six independent Gaussian distributed variates with shifted mean jy; and reduced variance (1 — 7 2 ), i.e. 
m i = 72/i + \/l — 7 2 /i- 

Equations (|23|l , ()24|) and (|25|l suggest a new excursion-set-based algorithm, in order to obtain the constrained eigenvalues 
of the matrix having components T a \H a (see Appendix [A] for the definition of the index a) - supplemented by the condition 
of a positive curvature H for the density field (or negative curvature, for voids). The procedure can be summarized as follows: 

(i) Draw six independent zero-mean unit-variance Gaussian distributed variates yi (i = 1, 6), and determine the components 
H a of the density Hessian matrix via Equation (|23|) . Compute the value of <th using (|A6|) . 

(ii) Calculate the eigenvalues f 1,^2, £3 of the previous Hessian matrix, and check if they are all positive (negative). If so, 
proceed to the next step, otherwise try again. This will guarantee the Hessian to be positive (negative) definite (i.e. the 
Hessian is a real-symmetric matrix), which is the condition for maxima (minima) of the field. Note that this step is clearly 
not required if we want to sample only p(T|H,7). 
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(iii) Draw other six independent Gaussian distributed variates k (i — 1,6), with mean zero and variance one, and determine 
the T a \H a components via Equation (|25[) . using the previously accepted values H a from (|23[) - while ot is determined via 
(|A6[) . Since we require the density Hessian to be positive definite, this means that we are effectively sampling the conditional 
probability p(T|H > 0, 7) - or p(T|H < 0, 7) for a negative definite Hessian. 

(iv) Calculate and store the constrained eigenvalues Q of the matrix having components T a \H a > 0. 

We leave in Appendix [B] some more insights on the main conditional formula (fT]), which readily follow from the previous 
Equations (|23[) . (|24[) and (|25[) . In the next sections, we will show a few applications of the new algorithm - with particular 
emphasis on the conditional distributions of shape parameters in triaxial models of structure formation (i.e. ellipticity and 
prolateness). Additional applications will be presented in forthcoming publications. 



4 CONDITIONAL DISTRIBUTIONS AND PROBABILITIES: NUMERICAL TESTS 

The algorithm illustrated in the previous section allows one to test and confirm several formulas derived by Rossi (2012), 
regarding conditional distributions and probabilities subjected to the extremum constraint. We present here results of the 
comparison (i.e. theory versus numerical), where for simplicity we set or = oh = 1; this corresponds to adopt 'reduced 
variables', as done in Rossi (2012). Note that we show explicitly the dependencies on a values in the following expressions, 
although we set them to unity in the mock tests. 



4.1 Individual conditional distributions 

The individual conditional distributions of shear eigenvalues, for a given correlation strength 7 with the density field (which 
encapsulate the peak/dip constraint), are given by (Rossi 2012): 

\/5 r 20 Ci / 9 Cf\ \/27r / 5 Ci ' 



p(Ci|7) = i2^L(ww CX H' w^fwJ a-7 2 ) 3 / 2 ex H~ w^rwJ (26) 

«*( - ""7==") [d - 7 2 ) - 204l + ^= «P( - 777^4 W - 

V y/l -72 a T J 17 ,J a^l ^1_ 7 2 H V 4(l-7 2 )cr2y V 2^/1 - 7 2 ctt/J 



The previous formulas show explicitly that the conditional distributions of shear eigenvalues are Doroshkevich-like ex- 
pressions, with shifted mean and reduced variance. Different panels in Figure [1] display plots of these distributions, contrasted 
with results from the algorithm presented in Section [3l where 500,000 mock realizations are considered. In particular, note 
the symmetry between p(£i|7) and p(Cs|7)- We expect the mean values and variances of those distributions to be: 

<Ci|7> = ^=-T\^f, 13 3 ~ 27 4(l- 7 2 ) (29) 

<C2|7> = 0, 4h^4(l-7 2 ) (30) 

2 13tt - 27 



<Ca|7> = — <Ci|T> = attF^tV^ -7 2 , <4,It = °Ci|-r = — ™ ct t(!-7 2 )- (31) 

V107T OU7T 

In the various panels, both theoretical and numerical expectations are reported and found to be in excellent agreement. For 
example, when 7 = 0.50, (Ci|7)num = 0.4644 while the expected theoretical value is (Ci|7)th = 0.4635, and for the variances 
""CiTt™ = 0-H02 where the theoretical expectation is o" 2 ^ = 0.1101. 

Obviously, to see explicitly the effect of the peak constraint, one needs to make cuts in £i from the partial conditional 
distributions Ai|£i, since = Aj — rj^i- In particular, expressions for p(Ai|H > 0, 7) can be derived from (|26H28[) via substitution 
of variables and integration over H > 0. For example, it is straightforward to obtain an analytic formula for p(As|H > 0,7). 
In Section [5] we will return on these issues in more detail, in connection with the conditional ellipticity and prolateness for 
dark matter halos. 

Another important distribution is p(A T | H |C3 > 0,7), with A T | H = K\ = (1 + (2 + (3 being the sum of the constrained 
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Figure 1. Individual conditional probabilities p(Ci I T) [top], p(C2\l) [middle], and p(C:s|7) [bottom] of the initial shear field in the 
peak/dip picture. Solid curves are Equations H26I I, 12711 , 128p . while histograms are obtained from 500,000 mock realizations via the 
algorithm described in Section \3\ Each panel displays the numerical and theoretical expectations for the mean values and variances of 
these distributions; their agreement is excellent. 



eigenvalues, namely the probability distribution of A T | H confined in the regions with £3 > (i.e. all positive eigenvalues): 

1/ ^ n \ 75%/5 A T |h / 9 Arjg A , , qo ^ 

p(A T|H |C3>0 )7 ) = - — ^— exp(- I g , (1 _ 72) )+ (32) 

+ I cy J A ^iH )\ eli ( ^ A t|h \ crf / VTOAte N1 

4V2^ctt\/1-7 2 V 2c4(l- 7 2)n \4cttx/1-7 2 ^ \2o-iVl - 7 2 '- 1 

This distribution is expected to have mean value 



(A T | H |C3 > 0,7) = ^12(3^6-2) a-rVl-7 2 =i 1.66 aryT^f. (33) 
144 y/ir 

It is also easy to see that the maximum of p(A T | H |C3 > 0,7) is reached when A T | H — 1.5 ot\/1 — 7 2 , and with a more 
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-4 -2 2 4 -4 -2 2 4 -4 -2 2 4 -4 -2 2 4 

(A T|H |C 3 >0,7) 

Figure 2. Conditional probability distribution p( A T | jj IC3 > 0)7)- Solid lines are obtained from Equation (132 b . for different values of 
the correlation parameter 7, while histograms are drawn from 500,000 realizations via the algorithm described in Section [3] Note the 
excellent agreement between numerical results and theoretical predictions for the corresponding mean and rms values, reported in the 
various panels. In addition, solid vertical lines show the expected maximum value of each distribution, while dotted lines represent their 
corresponding mean value. 



elaborate calculation its variance can be estimated analytically. The result is: 

<& TIH |Ca>o,7 = 1^(1 ~ 7 2 ) [arctan(v / 5/2) + arctan(^5) - ^ - ^(3^ - 2) 2 ] ~ 0.31 <4(1 - 7 2 ). (34) 

Figure [2] confirms the previous relations, by contrasting Equations (|32[) . (|33[) and (|34l) with numerical results from 500,000 
realizations from the algorithm presented in Section [3] We find good agreement, and recover correctly the expected mean 
values and variances of the distributions. For example, when 7 = 0.25, we find (A T | H |C3 > 0, 7 } n um = 1.6090 while the 
expected theoretical value is (A T | H |^3 > 0, 7 } t h = 1.6040, and for the rms values we find o"a™ h |C3>0 7 ~ 0.5314 where the 
expected theoretical value is o"a T | H |c 3 >0 7 — 0.5399. In the absence of correlations between the potential and density fields 
(i.e. when 7 = 0), all the previous expressions reduce consistently to the unconditional limit of Lee & Shandarin (1998). 



4.2 Distribution of peak heights 

An important quantity which plays a major role in peaks theory (Bardeen et al. 1986) is the distribution of peak heights, 
namely p(f|H > 0,7), where v — St/<jt and the overdensity St is the trace of the shear tensor - defined in Section [21 
One can easily obtain a simplified analytic expression for this distribution as follows. First, normalize Equation (|12l) by <jt, 
defining N = At|h/o"t and x = <5h/o"h, so that N = v — 73;. In this framework, x is the peak curvature of Bardeen et 
al. (1986) if Gaussian filters are used for computing the moments of the smoothed power spectrum (Equation IA6|) . Since 
p(N\j) is simply a Gaussian with mean zero and variance (1 — ■y 2 ) - see Rossi (2012) - and because of Equation (|12[) . then 
clearly p(y\x,j) is also Gaussian, with mean 71 and reduced variance (1 — 7 2 ). Note in fact that p(y, x^) is a bivariate 
Gaussian, because both v and x are normally distributed random variates. This implies that (v\x, 7) = 7s, and clearly 
(z/|H > 0,7) = 7(z|£3 > 0) = 7(^|A 3 > 0) (Bardeen et al. 1986) - where (i/|A 3 > 0) = 1.6566 according to Equation (21) of 
Lee & Shandarin (1998). Starting from a bivariate Gaussian with the previous expected mean value, it is then direct to derive 
a fairly good approximation for the distribution of peak heights, namely 

PMH > o, 7 , . J= « P [<^i£] [. + m 

where 

X=^[f^(3^-2)-l] 7 ^ 0.867. (36) 

In essence, p(f|H > 0, 7) is a Gaussian with shifted mean modulated by the role of the peak curvature, and it consistently 
reduces to a zero-mean unit-variance Gaussian distribution when 7 = - as expected. Additionally, the previous distribution 
can be shown to be equivalent to p(f)[/ °°p(^3|7, v)&(,3 / p{(,3 > 0)], where the term in square brackets quantifies the effect of 
the peak constrain - and so the scale at which the difference between peaks and random positions would appear. 

Figure [3] shows plots of p(^|H > 0,7) for different values of 7. In the various panels, vertical lines are the expected 
mean values; we find good agreement between numerical estimates and theoretical expectations. For example, when 7 = 0.25, 
(f|H > 0, 7) aum = 0.4146 while (f|H > 0, 7) t h = 0.4142. Solid curves in the figure are drawn from Equation (f35j) . an 
approximation which is particularly good for lower values of the correlation strength. 
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Figure 3. Distribution of peak heights, p(v\H > 0,7), for different values of the correlation parameter 7 - as specified in the various 
panels. The numerical values for the mean of the distribution as a function of 7 are indicated in each panel, and also marked with vertical 
lines; they are in good agreement with theoretical expectations. Solid curves are drawn from Equation (135 D - 



5 CONDITIONAL ELLIPTICITY AND PROLATENESS 

The theoretical framework outlined in Section[2] along with the algebra leading to the new computational procedure presented 
in Section [3l provides a natural way to include the peak/dip constraint into the initial shape distributions of halos and voids. 
This is a key aspect for implementing an excursion set algorithm for peaks and dips in Gaussian random fields (which will be 
discussed in Section [6}, and has direct applications in the context of triaxial models of structure formation. In particular, our 
description combines the formalism of Bardeen et al. (1986) - based on the density field - with that of Bond & Myers (1996) 
- focused on the shear field. 

In what follows, first we introduce some useful definitions regarding shape parameters (ellipticity and prolateness); we 
then provide a new expression for the joint conditional distribution of halo shapes given the density peak constraint, derived 
from Equation (jTTj) , which generalizes some previous literature work. The formula has important implications for the modeling 
of non-spherical dark matter halos, in relation to their initial shape distribution. Finally, we briefly discuss how the shear 
ellipticity and prolateness (and so the initial halo shapes) are modified by the inclusion of the density peak condition. 

Our analysis is based on the constrained eigenvalues of the initial shear field (Equations I11H 91. and we use the new 
algorithm to support our description. While here we only consider the case for dark matter halos, it is straightforward to 
deal with voids. In a forthcoming companion publication, we will present a more detailed study on the morphology of both 
halos and voids in the peak/dip picture, where we also investigate the modifications induced by primordial non-Gaussianity 
on their shapes (see also Section [7]) . 



5.1 Shape parameters: general definitions 

In triaxial models of halo formation, such as the ellipsoidal collapse (Icke 1973; White & Silk 1979; Barrow & Silk 198; Kuhlman 
et al. 1996; Bond & Myers 1996), it is customary to characterize the shape of a region by its ellipticity and prolateness. In 
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particular, in the 'peak-patch' approach of Bond & Myers (1996), the shape parameters are associated with the eigenvalues 
of the shear field ( Ai in our notation) . Considering the mapping ( Ai , A2 , A3 — > eT , Pt , St ) , where the subscript T refers to the 
tidal tensor, we can write: 

r _ , \ 1 \ 1 \ „ Ai — A 3 Ai + A 3 — 2A 2 _ A 2 — A 3 

o T = fci = Ai + A 2 + A3, e T = — 77 , Pt = ^7 = e T = , (37) 

zot 2ot ot 

where eT and pt are the 'unconditional' shear ellipticity and prolateness. In particular, if St > then eT > and therefore 
— e T < Pt < e T . 

In the 'peaks theory' of Bardeen et al. (1986), instead, ellipticity (en) and prolateness (pn) are associated with the 
eigenvalues of the density field (£i in our notation). Hence, given the mapping (£1, £2, £3 — > en , prr , <5h) where now the subscript 
H refers to the Hessian of the density field, one has: 

r l £l — £3 £l + £3 — 2^2 £2 — £3 . . 

Sh = hi = £1 +£2 + £3, e H = - ^ , p H = ^ — s e H - (38) 

2oh 2Ah oh 

As in the previous case, if Sh > then eH > and — eH ^ ph ^ eH- 

It is therefore natural to consider the 'conditional' mapping (Ci, C21 C3 ~~ * ^tihi -Ft|Hj A T | H ) an d define: 

A T | H = Ki = Qi + £2 + Qs, ^t|H = 77a ' "t|h = ^rr = ^t|h t , (39) 

ZL±t\-R Z£±t\h ^T|H 

where £i are the constrained shear eigenvalues - given by Equation (19). In what follows, we will refer to £tih an d Pt\h 
as to the conditional shear ellipticity and prolateness, respectively; we will also relate these quantities to the unconditional 
expressions previously introduced, involving (bt,Pt) and (eH , Ph) • 



5.2 Joint conditional distribution of shear ellipticity and prolateness in the peak/dip picture 

Combining ()37|) with Doroshkevich's formula (|20|) . it is straightforward to derive the 'unconditional' distribution of eT and 
Pt given St, p(ct,Pt|5t), related to the gravitational potential. In particular, 

9 (St , eT, pt) = p(oT)ff(eT,PT|f>T) (40) 

where p(o~t) is simply a Gaussian with zero mean and variance <r T , and 

15 3 / St \ 5 2 2 T 5 / St \ 2 2 2I 
g(e T ,PT^T) = - 7== — e T (e T - p T ) exp - - — (3e T +p T ) . (41) 

3V 107T V <TT / L 2 \ OT ' J 

Equation (|40[) implies that the unconditional joint distribution of shear ellipticity and prolateness is independent of that of 
the overdensity St- Similarly, inserting Q38p in Doroshkevich's formula (|2ip . one obtains 

g(Sn,eH,Pn) = p(8 K )g(e K , p H |fe) (42) 
for the quantities related to the density Hessian, where p(&) is a Gaussian with zero mean and variance u^, while 

g(e H ,PH[fa) = o 1 E— (— ) e H (e|-p|) exp[- |f— ) (3e|+p|)l. (43) 

Along the same lines, combining (|39|) with (|1 If) and using the definitions introduced in the previous section, it is direct 
to obtain 

G(A T | H , £ T | H , P T |h|7) = p(A T | H |7)G(£ T |h, P T |h| A T | H , 7) (44) 
where p(A T | H |7) is a Gaussian distribution with mean zero and variance cr T (l — 7 2 ) - i.e. Equation (57) in Rossi (2012) - and 

G(g T| H,PT|HlAT|H,7) = ^(^ ^ (45) 

The previous expression (|45[) is the joint conditional distribution of shear ellipticity and prolateness, given A T | H and a 
correlation strength 7 with the density field. Clearly, G(Et\h, Pt\h\ At|h, 7) is also independent of the distribution of (A T |h|7)- 
Equation (|45|) generalizes the standard joint distribution of halo shape parameters to include the peak constraint, and is another 
main result of this paper. Note that, in the absence of correlation (i.e. when 7 = 0), Et\h — > eT, Pt|h ~> Pt, A T | H — > St, so 
that ()45p reduces consistently to the 'unconditional' Doroshkevich's limit (|4ip . 

From this joint conditional probability function, it is possible to derive the partial conditional distributions for the 
shear ellipticity and prolateness given the peak constrain (which we will discuss next), and their expressions at density peak 
locations, when the condition H > is satisfied. 
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Figure 4. [Top] Distribution of (£1 — C3 It) ; f° r different values of the correlation strength 7 — as indicated in the panels. Its shape is 
directly related to the shear conditional ellipticity in the peak/dip picture - see Equation II46I I, [Bottom] Distribution of (£2 — C3 It)> f° r 
different values of the correlation strength 7. This distribution is instead related to the shear conditional prolateness in the peak/dip 
picture - see Equation l|47| l. The numerical values for the mean and variance of the distributions, as a function of 7, are also indicated 
in the figures. In both cases, solid curves are obtained via numerical integrations starting from the joint conditional distribution of 
eigenvalues {TTJ . Histograms are drawn from 500,000 realizations, using the algorithm described in Section \3\ 



5.3 Initial distributions of triaxial halo shapes at density peak locations 

We can readily express the conditional distributions of shear ellipticity and prolateness in the peak/dip picture in terms of 
the unconditional quantities (|37|) and ()38[) . This is done simply by recalling that the constrained shear eigenvalues are given 
by Q — Ai — r;£i, according to Equation (19). 
It is then direct to obtain: 

t? a (Ci ~ Cs) _ (Ai - A3) (£1 - £3) x x ,.„■. 

-Et|hA T |h = ^ = 2 77 2 = ~~ f?o H e H . (46) 

The previous relation implies that, for a given A T | H , the conditional shear ellipticity will have its mean value shifted by 
the presence of the peak constraint; the entity of the shift has to be ascribed to the additional factor r/Snen, which is given 
by the density ellipticity (essentially, the latter term quantifies the role of the peak curvature). The top panel of Figure [4] 
shows the distribution of (£1 — C3I7) f° r different values of the correlation strength 7, namely the combination of constrained 
shear eigenvalues which controls the conditional ellipticity - according to Equation (|46[) . Histograms are drawn from 500,000 
realizations using the algorithm described in Section [3] while solid curves are obtained via numerical integration - starting 
from the joint conditional distribution of eigenvalues pip . In the various panels, we also provide the numerical values for the 
mean and variance of the distribution, as a function of 7. 
Similarly, 

pA — (Ci + Cs - 2C 2 ) _ (A! + A 3 -2A 2 ) (ft+fr-2fr) _ 

i-T| H A T | H — - — 77 - — 6tPt - ^OHPH = ii T |H^T|H — (C2 — (,3)- (47) 

This expression shows that the conditional shear prolateness can be obtained from the conditional shear ellipticity (|46|l and 
the combination (£2 — C3) °f the constrained shear eigenvalues. The bottom panel of Figure [4] displays the distribution of 
(C2 — C3I7); f° r different values of the correlation strength 7. Again, solid curves are derived from a numerical integration of 
Equation (fTT) . 

Of particular interest are also the distributions <?(eT|H > 0,7) and g(pT\U > 0,7), which will provide the initial triaxial 
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Figure 5. Distributions p(Ci — C'i\C'A > 0, 7) [top] and p(C2 — C3IC3 > 0,7) [bottom] for different values of the correlation strength 7. In 
both figures, histograms are drawn from 500,000 realizations using the algorithm described in Section [3] while solid curves are obtained 
via numerical integrations of the joint conditional distribution of eigenvalues lilt . In the various panels, we also provide the numerical 
values for the mean and variance of the distributions, as a function of 7. In particular, their mean values at 7 = allow one to quantify 
the shift in the mean of the shear ellipticity and prolateness caused by the peak constraint. 



shapes of dark matter halos at peak locations. They are related to the previous expressions, and can be readily obtained 
within the outlined formalism; clearly, they will reduce to the standard (or unconditional) shear ellipticity and prolateness 
when 7 = 0, in the absence of correlation. We will present a dedicated study focused on the morphology of halos and voids 
in the peak/dip picture, where we characterize in detail those distributions, especially in relation to the ellipsoidal collapse 
model. We anticipate here that several analytic and insightful results on their shapes can be derived, starting from (|45[) and 
using (|46[) an d <|47[) . For example, according to (|46[) . when H > (which is equivalent to impose the condition £3 > 

on an ordered set of density Hessian eigenvalues), then the quantity (£1 — £3 1^3 > 0} will essentially provide by how much the 
mean value of the shear ellipticity has been shifted by the peak constrain. For the prolateness the situation is slightly more 
complicated, but according to (|47|) the shift can be derived by the additional knowledge of (£2 — §3 1 ^3 > 0). To this end, Figure 
Oshows the distributions of (£1 — C3K3 > 0,7) an d (C2 — C3IC3 > 0,7), for different values of the correlation strength 7 (note 
that the quantities just discussed are their mean values when 7 = 0). As in the previous plot, histograms are drawn from 
500,000 realizations using the algorithm described in Section [3] while solid curves are obtained via numerical integrations - 
starting from the joint conditional distribution of eigenvalues (|11|) . In the various panels, we also provide the numerical values 
for the mean and variance of the distribution, as a function of 7. 

As mentioned before, we will return in more detail on these distributions in a companion publication, to quantify the 
impact of the peak constrain on their shapes. We will also make the connection with the work of Bardeen et al. (1986) more 
explicit - see their Appendix C. It will be also interesting to include the role of the peculiar gravity field in our description, 
along the lines of van de Weygaert & Bertschinger (1996), as well as to extend the work of Desjacques (2008) on the joint 
statistics of the shear tensor and on the dynamical aspect of the environmental dependence, within this formalism. 



6 EXCURSION SET APPROACH FOR PEAKS AND DIPS IN GAUSSIAN RANDOM FIELDS 

The ability to distinguish between random positions and peaks/dips, contained in Equations (JXJ) or fTTJ and achieved by the 
algorithm presented in Section [3] is the key to implement an excursion set model for peaks and dips in Gaussian random 
fields. Indeed, the primary motivation (and one of the main strengths) of the proposed sampling procedure resides in its 
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direct inclusion into the excursion set framework (Epstein 1983; Peacock & Heavens 1990; Bond et al. 1991; Lacey & Cole 
1993). In essence, because in our formalism the peak overdensity is simply the trace of the conditional shear tensor (recall 
the definitions in Section [2j, and since in triaxial models of collapse the initial shape parameters are just combinations of 
the shear eigenvalues (Bond & Myers 1996), our prescription provides a direct way to generate the distribution of initial 
overdensities under the conditions that they are peaks/dips (i.e. the distribution of peak heights, when H > 0) - along 
with the corresponding conditional distribution of initial shapes (see Section [S]). In this respect, it is then straightforward to 
include this part in standard excursion set algorithms, as those used for example in Chiueh & Lee (2001), Sheth & Tormen 
(2002) or Sandvik et al. (2007) to compute the mass function, or in Rossi, Sheth & Tormen (2011) to describe halo shapes; 
the 'peak/dip excursion-set-based' algorithm is also useful for computing the halo bias assuming triaxial models of structure 
formation (i.e. ellipsoidal collapse). The only main conceptual difference is the pre-selection of peak/dip locations, instead 
of random positions in the field. The mass scale of the peak will then be fixed by finding the proper ctt which satisfies the 
combination of (5t, e-r,PT|H > 0,7), assuming some structure formation models - as for example the ellipsoidal collapse. 

This is particularly useful because the excursion set theory is a powerful tool for understanding various aspects of the full 
dynamical complexity of halo formation. Perturbations are assumed to evolve stochastically with the smoothing scale, and the 
problem of computing the probability of halo formation is mapped into the classical first-passage time problem in the presence 
of a barrier. A very elegant reformulation of this theory has been recently proposed by Maggiore & Riotto (2010a,b,c), who 
made several technical and conceptual improvements (i.e. no ad hoc absorbing barrier boundary conditions, account for non- 
markovianity induced by filtering, unambiguous mass association to a smoothed scale, etc.) by deriving the original excursion 
set theory from a path integral formulation - following a microscopical approach. These authors also noted that the failure 
of the standard excursion set approach may be related to the inadequacy of the oversimplified physical model adopted for 
halo formation (either spherical or ellipsoidal) , and propose to treat the critical threshold for collapse as a stochastic variable 
which better captures some of the dynamical complexity of the halo formation phenomenon. Even so, they find that the 
non-markovian contributions do not alleviate the discrepancy between excursion set predictions and iV-body simulations. 

In addition to the problems pointed out by Maggiore & Riotto (2010a, b,c), there is also the fact that - in its standard 
formulation - the excursion set approach is unable to differentiate between peaks/dips and random locations in space - i.e. all 
points are treated equally. However, since local extrema are plausible sites for the formation of nonlinear structures and there 
is a good correspondence between peaks in the initial conditions and halos at late times, it may be important to differentiate 
between those special positions in space. The algorithm proposed here goes in this direction, since it allows one to pre-select 
those special points in space around which halos form (peaks), and not just random locations - and permits to associate 
their corresponding initial shape distribution: in essence, the computational procedure selects a special subset, among all the 
possible random walks considered in the standard excursion set procedure. 

Therefore, the 'peak/dip excursion-based' algorithm can be used to study the mass function of halos and their triaxial 
shapes at peak/dip positions, and also the halo bias. In fact, our prescription allows one to generate the initial distributions 
of overdensity, ellipticity and prolateness (i.e. shape parameters) at a scale set by the variance or, with the constraint that 
St is a peak (i.e. the condition H > on the Hessian). One can then just evolve this initial conditional shape distribution 
by solving a dynamical equation of collapse, and study the final shape distribution (as done for example in Rossi, Sheth 
& Tormen 2011) but now at peak/dip locations. Alternatively, one can adopt a 'crossing threshold', since in the excursion 
set approach an halo is formed when the smoothed density perturbation reaches a critical value for the first time, and the 
problem is reduced to a first-passage problem in the presence of a barrier (i.e. if the overdensity exceeds a critical value, the 
random walk stops at this scale; if not, the walk continues to smaller scales). Moreover, our numerical technique can be easily 
integrated in Montecarlo realizations of the trajectories obtained from a Langevin equation with colored noise (i.e. Bond et 
al. 1991; Robertson et al. 2009) at peak/dip locations, and even for situations where the walks are correlated - in presence of 
non-markovian effects, along the lines of De Simone et al. (2011a,b). In this respect, our technique is general, since any kind 
of filter function can be readily implemented. All these lines of research are ongoing efforts, and will be presented in several 
forthcoming publications. 



7 CONCLUSIONS 

From the joint conditional probability distribution of an ordered set of shear eigenvalues in the peak/dip picture (Rossi 
2012; Section [2] Equations [T1 and 1 1 1 [) . we derived a new algorithm to sample the constrained eigenvalues of the initial shear 
field associated with Gaussian statistics at peak/dips positions in the correlated density field. The algorithm, described in 
Section[3] was then used to test and confirm several formulas presented in Rossi (2012) regarding conditional distributions and 
probabilities, subjected to the extremum constraint (Section |4)|. We found excellent agreement between numerical results and 
theoretical predictions (Figures [1] and [2}. In addition, we also showed how the standard distributions of shear ellipticity and 
prolateness in triaxial models of structure formation are modified by the constraint ( Section Figures [3] and [SJ), and provided 
a new expression for the conditional distribution of shape parameters given the density peak requirement (Equation I45[) . 
which generalizes some previous literature work. The formula has important implications for the modeling of non-spherical 
dark matter halo shapes, in relation to their initial shape distribution, and is directly applicable to the ellipsoidal collapse 
model (Icke 1973; White & Silk 1979; Barrow & Silk 1981; Kuhlman et al. 1996; Bond & Myers 1996). In particular, our novel 
description is able to combine consistently, for the first time, the formalism of Bardeen et al. (1986) - based on the density 
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field - with that of Bond & Myers (1996) - based on the shear field. Along the way, we also discussed the distribution of peak 
heights (see Figure [3]), which plays a major role in peaks theory (Bardeen et al. 1986). 

While the primary motivation of this paper was to illustrate the new 'peak/dip excursion-set-based' algorithm, and to 
show a few applications focused on the morphology of the cosmic web (following up, and complementing with some more 
insights, the theoretical work presented in Rossi 2012), the other goal was to describe how the new sampling procedure 
naturally integrates into the standard excursion set framework (Epstein 1983; Peacock & Heavens 1990; Bond et al. 1991; 
Lacey & Cole 1993) - potentially solving some of its well-known problems. In particular, in Section H5] we argued that the 
ability to distinguish between random positions and peaks/dips, encoded in the algorithm and in Equations (fTJ) and fTTJ 
derived from first principles, is indeed the key to implement a generalized excursion set model for peaks and dips in Gaussian 
random fields. This is the actual strength of the proposed computational procedure, since part of the failure of the original 
excursion set theory may be attributed to its inability to differentiate between random positions and special points (peaks) 
in space around which halos may form. 

To this end, our simple prescription can be used to study the halo mass function, halo/void shapes and bias at peak/dip 
density locations, in conjunction with triaxial models of structure formation. All these research lines are ongoing efforts, 
subjects of several forthcoming publications. The essential part is the characterization of the distribution of peak heights, and 
of the initial shape distribution at peak/dip locations f Equations [Tl [TT1 and 145 [1. 

The algorithm presented in this paper offers also a much broader spectrum of applications. This is because, as pointed 
out by Rossi (2012), the fact that the eigenvalues of the Hessian matrix can be used to discriminate between different types of 
structures in a particle distribution is fundamental to a number of structure- finding algorithms, shape-finders algorithms, and 
structure reconstruction on the basis of tessellations. For example, it can be used for studying the dynamics and morphology of 
cosmic voids - see for example van de Weygaert & Platen (2011), Bos et al. (2012), Pan et al. (2012), and the Monge-Ampere- 
Kantorovitch reconstruction procedure by Lavaux & Wandelt (2010) - and in several observationally-oriented applications, 
in particular for developing algorithms to find and classify structures in the cosmic web or in relation to its skeleton (Sahni 
et al. 1998; Schaap & van de Weygaert 2000; Novikov et al. 2006; Hahn et al. 2007; Romano-Diaz & van de Weygaert 
2007; Forero-Romero et al. 2009; Zhang et al. 2009; Lee & Springel 2010; Aragon-Calvo et al. 2010a,b; Platen et al. 2011; 
Cautun et al. 2012; Hidding et al. 2012). Another application is related to the work of Bond, Strauss & Cen (2010), who 
presented an algorithm that uses the eigenvectors of the Hessian matrix of the smoothed galaxy distribution to identify 
individual filamentary structures. In addition, since galaxy clusters are related to primordial density peaks, and there is a 
correspondence between structures in the evolved density field and local properties of the linear tidal shear, our theoretical 
framework provides a direct way to relate initial conditions and observables from galaxy clusters. 

Other intriguing connections involve topological studies of the cosmic web, the genus statistics and Minkowski functionals 
(Gott et al. 1986, 1989; Park et al. 1991, 2005; Matsubara 2010), and the possibility to address open questions regarding the 
origin of angular momentum and halo spin within this framework; this is because the dependence of the spin alignment on the 
morphology of the large-scale mass distribution is due to the difference in shape of the tidal fields in different environments, 
and most of the halo properties depend significantly on environment, and in particular on the tidal field - i.e. the environmental 
dependence of halo assembly time and unbound substructure fraction has its origin from the tidal field (Wang et al. 2011). 
It will be also interesting to explore how the new formalism proposed here can be used to study halo spin, shape and the 
orbital angular momentum of subhaloes relative to the LSS, in the context of the eigenvectors of the velocity shear tensor 
(see the recent study by Libeskind et al. 2013). In addition, the more complex question of the local expected density field 
alignment /orientation distribution as a function of the local field value (Bond 1987; Lee & Pen 2002; Porciani et al. 2002; 
Lee, Hahn & Porciani 2009; Lee 2011) can be addressed within this framework, and is the subject of future studies. 

On the theoretical side, we note that our algorithm is restricted to one scale (i.e. peaks and dips in the density field, as in 
Bardeen et al. 1986), but the extension to a multiscale peak-patch approach along the lines of Bond & Myers (1996) is doable 
and subject of ongoing work. As argued in Rossi (2012), this will allow to account for the role of the peculiar gravity field 
itself, an important aspect not considered in our formalism but discussed in detail in van de Weygaert & Bertschinger (1996). 
Including all these effects in our framework and translating them into a more elaborated algorithm is ongoing effort, and will 
allow us to make the connection with the multiscale analysis of the Hessian matrix of the density field by van de Weygaert & 
Bertschinger (1996) and Aragon-Calvo et al. (2007; 2010a,b). It will also allow us to incorporate the distortion effect of the 
peculiar gravity field in our initial distribution of halo/void shapes (Section[SJ). The natural extension of the peak/dip picture 
for the initial shear to non-Gaussian fields is also ongoing effort, along with some other broader applications in the context of 
the excursion set model - for example in relation to the hot and cold spots in the Cosmic Microwave Background, including 
the effects of /NL-type non-Gaussianity on their shapes (i.e. Rossi, Chingangbam & Park 2011) - which will be presented in 
forthcoming publications. 
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APPENDIX A: ESSENTIAL NOTATION 

We summarize here the basic notation adopted in the paper, which is essentially the same as the one introduced by Rossi 
(2012) - with a few minor changes to make the connection with previous literature more explicit. In particular, let denote 
the displacement field, $ the potential of the displacement field (i.e. the gravitational field), F the source of the displacement 
field (i.e. the density field). Both F and $ are Gaussian random fields, the latter specified by the matter power spectrum 
P(k), with k denoting the wave number and W(k) the smoothing kernel. Use Ty for the shear tensor (its eigenvalues are 
Ai, A2,Aa), Hij for the Hessian matrix of the density F (with eigenvalues £1,^2, £3), Jij for the Jacobian of the displacement 
field, where i, j — 1,2, 3. Indicate with q the Lagrangian coordinate, with x the Eulerian coordinate, where x(q) = q + ^(q). 
Clearly, 

r f ^ 9xi . „ m <9 2 <E> Tr d 2 F _ , Aa*; ^3 2 $ 



The correlations between the shear and density Hessian are expressed by 

15 



2 

(TijTki) = -y(Mm + <5ik<$ji + <5n<5jk) (A2) 



15 

where a'^ = S2 = cr 2 , — Sc, = a 2 , Tth = S4 = ct 2 , 5ij is the Kronecker delta, and 



2 

(tfijtfki) = ^| (5,Ai + <5ik5ji + <5ii<5 jk ) (A3) 



(Hjtf kl ) = ^(Mki + <5ik<5ji + 5ii<5 jk ) (A4) 



s n = ^ I k" />(/,-! ir-'(//xi/,- i.v.vi 
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l 

2^2 



f GO 

/ fc 2(j+1) P(jfe) W 2 (k) dk = S; 
Jo 



2(j + l)' 



(A6) 



In the main text we prefer to use ot and oh, rather than the more familiar an and a 2 , for their intuitive meaning. In particular, 
the subscript T always indicates that a quantity is related to the shear field, while the subscript H denotes a quantity linked 
to the Hessian of the density field. 

The shear and density Hessian Tij and Zfij are real symmetric tensors, so they are specified by 6 components. 
Whenever necessary, we label those components with the symbols a or j3 to indicate the various couples, where a,/3 = 
(1, 1), (2, 2), (3, 3), (1, 2), (1, 3), (2, 3). It is also useful to introduce the vectors T and H, derived from the components of their 
corresponding tensors, i.e. T = (Tn, T22, T33, T12, T13, 223) and H = (-H11, H22, -H33, #12, #13, #23)- The constrained eigenval- 
ues of the matrix having components T a \H a will be indicated with £1, ("2, (3. As shown in Rossi (2012), the covariance matrix 
of the joint probability distribution of T and H is simply 

_t_ / erfA r TH A\ 
15 VrraA ctjjA J 

where 



y= ({T a T a ) (T a H p )\ 
\(HpT a ) (HpHp)) 



(A7) 



B 




(A8) 



with I a (3 x 3) identity matrix and 
'reduced' correlation: 



a (3 x 3) null matrix. Finally, an important 'spectral parameter' often used here is the 



7 = Fth/otOH = 



&0&2 



(A9) 



which plays a crucial role in peaks theory (i.e. Bardeen et al. 1986). If Gaussian filters are used in (|A6[) . then our 7 is the 
same as the one introduced in Bardeen et al. (1986) - specified by their Equation (4.6a). Note that in the main text we also 
define r\ — 7<tt/o"h; if one adopts reduced variables (i.e. T a and H a normalized by their corresponding rms values <7t and 
oh), clearly 77 = 7. 



APPENDIX B: INVARIANTS FROM THE CONDITIONAL FORMULAS 

The algebra presented in Section [3] allows one to gain more insights into the joint conditional distribution of eigenvalues in the 
peak/dip picture (Equation [l}. For simplicity, in what follows we consider 'reduced' variables, so that the various components 
of T and H are now normalized by their corresponding rms values (<tt and <th, respectively). With some abuse of notation, 
we omit the tilde symbol (used instead in Rossi 2012) to distinguish between normalized and unnormalized quantities. It is 
then possible to characterize and study the properties of the first few elementary symmetric functions of degree n for the 
density Hessian, the shear tensor, and the conditional shear tensor - along the lines of Weyl (1948), Doroshkevich (1970), 
Sheth & Tormen (2002) and Desjacques (2008). In particular, it is direct to note that, using Equation (|23[) and considering 
six independent Gaussian random variates j/i (i = 1, 6) represented by the six-dimensional vector y, one obtains that the first 
two classical invariants 

hi = Hu + H22 + H 33 = -yi (Bl) 

til = h\ - 3/12 = \{yl +yt + yl + y% + yl) (B2) 
5 

are independent. This fact implies that 

P(H) 4» e " S " /l5PW * ) ' (B3) 



hence p(H) is the product of two independent distributions, where in particular p(h\) is a one-dimensional Gaussian with mean 
zero and unity variance. Note also that p(H)dH = p(y)dy, where p(y) = Ilf = i 9i 18 simply the product of six independent 
one-dimensional zero mean unit variance Gaussians gi. 
Similarly, for the shear tensor T, one obtains that 

ki = Tn + T 22 + T 33 = -21 (B4) 

k\ = k\ - 3k 2 = -(22 + z\ + z'i + z\ + zl) (B5) 

are also independent, with zi (i = 1,6) other six Gaussian random variates represented by the six-dimensional vector z. 
Therefore, 

^ (T) ^^fc e " 5fci/2 ^ (&lMfc3) (B6) 
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with p(fci) a one-dimensional Gaussian distribution. Note again that p(T)dT = p(z)dz, where p(z) = Y\f =1 gi is the product 
of six independent one-dimensional zero mean unit variance Gaussians g\. 

Following the previous logic, one would naturally expect that the quantities K\ and K 2 = K\ — 3A2 , defined in the main 
text (see Equation [2j , should also be independent. Indeed, it is direct to obtain that (Ravi Sheth, private communication): 



K x = =-(mi-72/i) = - v / l-7 2 h (B7) 
K\ = K\ - 3K 2 

r[(m 2 - 7?/2) 2 + (m 3 - 7J/3) 2 + (m.4 - r yyi) 2 + (m 5 - 7y 5 ) 2 + (m 6 - 72/6) 2 ] 



5 L 

= (1 ~/ ] {ll + ll + li + ll + i§), (B8) 

where h (i — 1,6) are other six independent Gaussian distributed variates with mean zero and unity variance, while m\ 
(i — 1, 6) are six Gaussian distributed variates with shifted mean 73/i and reduced variance (1 — j 2 ) , i.e. rci\ — jyi + \/l — J 2 li- 
Hence, the joint conditional distribution of eigenvalues in the peak/dip picture (Equation [TJ can be written as the product of 
two independent distributions as follows: 

-A'J/[2(l- 7 2 )] -,.3 5if|/[2(l- 7 2 )] 

P(T|H,7) = VMi _ i2) 8VW7t5/2 (L _^ ^(A- 1 |7)P(^3|7). (B9) 

This latter expression clearly shows that the distribution of constrained A'i = At|h is independent of the distribution of 
the constrained angular momentum A"f, where in particular p(Ki\"f) is a Gaussian with zero mean and variance given by 
(1 — j 2 ), while p(Ks) is a chi-square distribution with five degrees of freedom. Once again, note that p(T|H, 7)d(T|H) = 
p(m|y, 7)d(mjy), where now 

JL „(mi-7!/i) 2 /2(l-7 2 ) JL 

p(m|y,7) = TT- = TT^j (BIO) 

namely, p(m|y,7) is now the product of six independent one-dimensional Gaussians U with shifted mean 71/i and reduced 
variance (1 — r y 2 ) - represented by the six-dimensional vector m|y. 



